Statistics of extinction and survival in Lotka-Volterra systems 
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We analyze purely competitive many-species Lotka-Volterra systems with random interaction 
matrices, focusing the attention on statistical properties of their asymptotic states. Generic features 
of the evolution are outlined from a semiquantitative analysis of the phase-space structure, and 
extensive numerical simulations are performed to study the statistics of the extinctions. We find 
that the number of surviving species depends strongly on the statistical properties of the interaction 
matrix, and that the probability of survival is weakly correlated to specific initial conditions. 
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I. INTRODUCTION 

Systems of interacting biological species evolve through the long, slow and intricate process of natural selection 
0. Usually, the result of this process is so complex that the dynamics of such webs of coevolving species can be 
successfully represented, within relatively short time scales, by means of a dynamical system with stochastic elements 
A standard mathematical model for the joint evolution of M biological species with spatially homogeneous 
densities ni{t) {i — 1,2, M) is the generalized Lotka-Volterra system H, 



ni{t) = ni{t) 



M 



(* = 1,2,...,M). (1) 



For large values of M , it is reasonable -as a phenomenological approach- to choose the parameters and K,ij at 
random from given probability distributions. Whithin this type of representation, the dynamics of coevolving species 
can be characterized by statistical properties over different realizations of parameter sets. 

There are two biological systems that can potentially involve a large number of coevolving populations. The 
first one is an ecological system in which each population corresponds to a different biological species, as usually 
interpreted in the theory of population dynamics [3jp|. The other situation is a system in which each population 
represents a genotype accessible to a given species pTT In this situation, the number of populations can be sensibly 
larger than in the case of interacting species. Although in both cases coevolution is presumably well described by 
Eqs. the probability distributions to be assigned to the random parameters Kij, which represent the interaction 
between populations, are not necessarily similar. In fact, in an ecological system of several coevolving species, mutual 
interactions can be of different types (competition, symbiosis, parasitisra) . Within a given species, instead, it is 
expected that the interaction is mainly competitive, as in logistic models g. 

It is well known that, in a system where many individuals compete for a resource, the dynamics leads to the 
extinction of some of them and to the survival of others. This is indeed a basic fact of evolution in the Darwinian sense. 
Though the generalized Lotka-Volterra model (|l|) has been explored in detail , it seems that a full characterization 
-either deterministic or statistical- of the conditions under which a population becomes extinguished or survives in 
the competition process has not been achieved. In this paper, we aim at analyzing this particular problem from a 
statistical viewpoint. 
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We consider a large number of coevolving species or genotypes, each of them consisting of a population of identical 
individuals with density ni{t). These populations are supposed to evolve according to the Lotka-Volterra model (|^), 
subject to purely competitive interactions, i.e. with k.^ > for any pair Since we aim at analyzing the statistical 
properties of the dynamics, these coefficients will be drawn at random from a given distribution and will remain 
quenched from the initial time. 

For simplicity, we take = 1 V i , indicating that in the absence of competition the dynamics of all the populations 
are identical. We are thus implicitly identifying these populations with the genotypes accessible to a given species. 
Within this condition -that it is not essential to our interest and could in fact be easily relaxed- Eqs.(|^) reduce to 



M 



1 



(i = l,2,...,M). 



(2) 



All the coefficients will be chosen at random from the same distribution p{k), such that p{k) = for k < 0. 

In the next section we outline the behavior of the dynamical system (^) in phase space, showing that its evolution 
proceeds along a series of "pseudo-extinctions", in which some of the densities ni(t) can attain very low levels during 
long periods but, eventually, they recover significative values. A threshold for these pseudo-extinctions -that become 
consequently true extinctions- is suggested by the biological context of the problem. This threshold is introduced in 
our numerical study of Eqs. (|^) in Section 3, where we focus the attention on the statistics of extinct and surviving 
genotypes and try to characterize their long-time behavior in terms of their inital conditions. Our results are discussed 
in Section 4. 



II. PHASE SPACE STRUCTURE 



The evolution of the dynamical system (||) can be described in terms of a semi-quantitative analysis of the corre- 
sponding phase-space topology, which is determined by the fixed points of (^) and the associated invariant manifolds. 
The equations for the fix-point coordinates n* read 



n*\=0 (i-l,2,...,M) (3) 



and have, generically, 2^ solutions. In fact, each solution to these equations can be characterized by the number 
M' of non-zero coordinates (M' = 0,...,M); let us call such a solution an M' -equilibrium. For a given choice of 
the coefficients Kij the number of different M'-equilibria is C(M, M') = M\/M'\{M - M')\. Therefore -disregarding 
pathological choices of Kij- the total number of fixed points is X^m' ^(-^i ■^') — "^^^ ■ 

Since n* stands for a density, meaningful equilibria among the 2^^ fixed points are those with non-negative coor- 
dinates. In Appendix A it is proven that, for random Kij, the probability that all the non-zero coordinates of an 
M'-equilibrium [M' ^ 0) are positive is 

P{M') = 2^-^^'. (4) 

We stress that it is essential to this result that nij > V i, j, i.e. that the system is purely competitive. For large M, 
the number of equilibrium points with non-negative coordinates will therefore be approximately given by 



ON A'/ 

-^^'C(A/,Af')«2(-j . (5) 



It is interesting to note that, if > Kmin V i, j, all the non- negative equilibria will be confined to a certain volume V 
in phase space, since n* < 1/Kmin. This volume shrinks rapidly for growing M, as F = K'^tl/M\, and the density 
of equilibrium points -most of which are situated on the surface of V , where some of the coordinates vanish- grows 
correspondingly. 

The stability properties of the fixed points of system (||) can be fully analyzed in some very special cases only. For 
instance, as could be expected for this logistic-like dynamical system, the 0-equilibrium [n* = V i) can be proven 
to be always unstable. Moreover, for a random choice of positive Kij, 1-equilibria are stable with probability equal to 
M~^. Finally, the Af -equilibrium (n* 7^ V i) is stable if is a symmetric matrix. 
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Though we cannot give a detailed characterization of the stabihty of all of the M'-equilibria, it can be argued that, 
for a random system and for large M' and Af, the eigenvalues of the linearized problem should follow a semicircular 
distribution |8|. A typical equilibrium point is thus linearly unstable and it has approximately the same number of 
positive and negative eigenvalues. Correspondingly, the number of unstable and stable invariant manifolds for each 
equilibrium is more or less the same. Since the mathematical structure of system (|^) prevents both the divergence of 
orbits and changes of sign in the densities rii [t) , the invariant manifolds of positive equilibria are necessarily bounded 
and mutually connected, defining homoclinic and heteroclinic orbits. Most of these orbits lie on the surface of the 
volume that contains the positive equilibria, where some of the densities are exactly equal to zero. 

In summary, the portion of the M-dimensional phase space of system (|^) meaningful to our problem is populated by 
a large set of fixed points -of the order of (3/2)*^ in number- most of them having positive and negative eigenvalues, 
i.e. being unstable. They are confined to a volume of order 1/M! and, typically, are found on the surface of such 
volume. These equilibria are highly interconnected through invariant manifolds which lie also on that surface and 
connect stable and unstable eigenvectors. The number of those manifolds should be of order M(3/2)^^. 

With these elements in hand, the evolution along a typical phase-space trajectory of the dynamical system (^) can 
be outlined as follows. From a generic initial condition, the orbit should soon approach one of the stable manifolds, 
and the system will be driven towards the corresponding equilibrium. It will spend some time in the vicinity of this 
equilibrium, but if this fixed point is not stable (i.e, if it has at least one unstable manifold, which -as we have argued- 
is the typical case) the orbit will finally leave that neighborhood, just to be drawn along one of the unstable manifolds 
of this first equilibrium point towards another equilibrium, that is expected to have in turn some stable and some 
unstable manifolds. The whole process will repeat itself and the system will wander in phase space, typically visiting 
the neighborhoods of a large number of unstable fixed points, until it eventually finds a stable equilibrium. This is 
reminiscent of the complex behavior of Boolean evolution models on random landscapes |^ , which -in contrast with 
Lotka-Volterra models- are however discrete (in space and time) and stochastic. 

We stress that, in wandering from one equilibrium to another, the orbit is expected to approach more and more the 
successive invariant manifolds that drive the dynamics of the system |l^. This implies, in particular, that the system 
will spend longer and longer periods in the immediate vicinity of those equilibria. Since, typically, the equilibria 
have some null coordinates, the corresponding densities will approach a vanishing state but — as the system escapes 
from each unstable equilibrium point — they can eventually recover appreciable values. As the evolution proceeds, 
a given density can therefore practically vanish during a rather long time, but can then increase and become again 
significative in the whole dynamics. We shall come back to these pseudo-extinctions in the next sections, to discuss 
their relevance in the numerical study of the system and its biological interpretation. 

Finally, it is worthwhile to remark that the existence of a stable equilibrium point, able to definitively attract an 
orbit, is in principle not guaranteed. Moreover, even if one or several stable points do exist, it is not insured that 
their bassins of attraction cover the whole space of initial conditions. The system could thus perform a chaotic orbit 
or become trapped in a limit cycle pO|JTl[ |. 

III. NUMERICAL ANALYSIS 

We have performed an extensive numerical investigation of system (^). Each realization consists of the numerical 
integration of the equations, after a random choice of the interaction matrix and the initial conditions. In all the 
results presented here, the interaction coefficients have been randomly chosen from a uniform distribution in the 
interval [kq ~ '^o + with Ak < kq, but we have tested that other probability distributions -always defined 
for K > 0- produce essentially the same results. Similarly, the initial densities have been uniformly distributed at 
random in the interval [0, nmax]. A proper rescaling of densities and time makes it possible to fix, without generality 
loss, kq = 1 and rimax = 1- The only parameter to vary in these distributions is therefore Ak. In the following we 
describe the dynamical behavior of (|2|) as drawn from our numerical calculations. 

A. Pseudo-extinctions and density threshold 

In Figure 1, we show the evolution of several typical densities in a system of M = 20 genotypes, for Ak = 0.5. Note 
that, to ease the appreciation of certain details, both the time and the density axis are logarithmic. In the inset the 
same curves are shown in a semilogarithmic plot, with logarithmic time axis. The phenomenon of pseudo-extintions is 
clearly seen in some of the curves. We have checked that, in some realizations, one or more densities can temporarily 
attain values as small as n '--^ 10^^^, and then grow to levels of the order of their initial values. The verification that 
pseudo-extinctions do occur -as predicted, in the previous section, from our analysis of the phase space structure- 
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points the attention to a further ingredient, which is not present in the model as described by (g), but has to be 
necessarily taken into account in a system where the variables are, actually, discrete. In fact, the population density 
of a species or a genotype confined to certain spatial domain 11 of volume Vq cannot be smaller than V^^, unless 
it vanishes. In a description in terms of densities, it is therefore necessary to fix a threshold [jl^ , below which the 
only value accessible to the density is effectively zero. Besides this biological argument for introducing a density 
threshold, we must stress that in our numerical calculations this element is also necessary to avoid spurious effects of 
finite computer precision on the results. A threshold no has been thus introduced as an additional parameter in the 
numerical calculations, in such a way that if a density attains a value lower than hq it is automatically set to zero. 

From the analytical viewpoint, it can be argued that the introduction of a threshold changes the stability of almost 
all the M'-equilibria. Roughly speaking, whereas without threshold an orbit could approach an equilibrium following 
a stable manifold just to leave it along an unstable one, with threshold the system can instead be "captured" by the 
equilibrium point if the orbit crosses the threshold. What was an unstable equilibrium becomes, in effective terms, a 
stable one. 

In order to illustrate the different behavior of systems with and without density threshold, we have chosen to analyze 
the evolution of the total density 

M 

4=1 

as a global characterization of the phase space dynamics. Figure 2 displays the evolution of N{t) for two systems of 
100 species. In both cases the interaction is defined by Ak = 1. One of the curves corresponds to the system without 
threshold (besides that imposed by the smallest representable number in the computer). The other one corresponds 
to the same system -with the same initial conditions and interaction matrix- with a threshold uq = 10^®. Note that 
the time scale is again logarithmic. It can be seen that, in this realization, the orbit of the first system follows the 
qualitative behavior described in the previous section. It passes near some equilibria -where N{t) remains practically 
constant- spending exponentially larger and larger times in their neighborhood. When a threshold is present, the 
system is always attracted to one of the new "stable" equilibria. In the case of Fig. 2, the system is captured at 
t w 100 by an equilibrium point that, although being unstable in the first case, acts now as a stable fixed point for 
this orbit. 



B. Statistics of survivals 

According to our simulations, the main feature in the long-time dynamics of system (||) -with or without threshold- 
is that it evolves to a situation in which most of the densities vanish -at finite times if no ^ or asymptotically if 
no = 0. This behavior can be identified with the extinction of the corresponding genotypes. In any case, for a large 
system, a variable number of surviving populations is found at long times. In Fig. 3, we show the distribution of 
the number of surviving populations for different values of Ak. Each curve was constructed from the results of 2000 
realizations in a system with M — 100 and no = 10~^. The final time in each realization was chosen so that a 
stationary state had been reached, which was checked to be a solution to Eqs. (||). 

The distribution of the number of surviving populations is generally a bell-shaped curve, its width and maximum 
depending on the probability distribution p{k). It can be seen in Fig. 3 that for An small enough the curve is 
relatively broad, and that wider interactions reduce the overall stability of the system, leading to a shift of the curve 
towards a situation where less species survive. The correlation between the maximum of the distribution of survivors 
and the width of p{k) is shown in the inset, in a log-log plot. Observe that for the smallest value, Ak = 0.02, the 
maximum of the distribution coincides with the total number of species in the system. 

It is worthwhile to note that the probability of having a non-negative M'-equilibrium, P{M') — 2^~*^ C{M, M') 
(cf. Section 2), is also a bell-shaped curve as a function of M' . This fact should however be considered only as 
indicative of the profile of the curves in Fig. 3. In fact, the probability that an initial condition approaches any 
M'-equilibrium depends not only on their number for a given M' but also on the size of their bassins of attraction, 
about which our phase-space analysis provides no information. 

C. Characterization of survivors 

A very natural question arises now, yet not easy to answer: which populations survive? Are they characterized 
by some particular initial condition, by some peculiarity in the interaction with the remaining populations? In the 
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study of a real ecological system, it would be desirable to give an answer to these questions in terms of quantities 
accessible from observations. It is thus reasonable to consider, since they are probably the easiest to measure, the 
initial densities and their initial time derivatives. Besides being accessible, these quantities characterize the initial 
interactive scenario: according to Eqs. (||), the density measures the effects of each population on itself, while its first 
time derivative accounts for the influence of the remaining species. 

We have found that an answer based on deterministic arguments cannot be given to such questions. According to 
the statistics collected from the simulations, we conclude that only a weak correlation exists between survival and the 
initial conditions. This correlation can be evidenced by calculating the distribution of final densities as a function of 
the initial one, providing thus a probabilistic answer to those questions. In Figs. 4 and 5 we show (full lines) the 
distribution of survivors as functions of initial values, for general asymmetric systems (k^ 7^ ^ji)- Fig. 4 shows that 
the probability of survival is almost uniform in the whole range of initial densities, with a slighty higher probability 
of survival for the largest ones. In Fig. 5 an associated distribution is shown: the number of survivors as a function of 
the initial derivative of the density. Here, an enhancement of the probability of surviving is seen around a relatively 
large (and negative) value of the initial time derivative. 

We have also found that this correlation between the final state and the initial condition is stronger in symmetric 
systems {nij = Hji). In Figs. 4 and 5, the same functions are shown as dashed lines for symmetric systems. In 
this case, there is an enhancement in the probability of survival for those species that start with a higher initial 
population. Although considering a purely symmetric interaction matrix is irrelevant from the biological point of 
view, these results suggest that the statistical correlation between initial and final states -and, in particular, between 
initial conditions and survival probability- can depend in a rather strong manner on additional constraints in the 
interaction matrix. 



IV. CONCLUSIONS 

We have analyzed a dynamical system that represents the evolution of many species coupled by Lotka-Volterra 
interactions. The study has been restricted to systems where the interaction is purely competitive. This dynamical 
system describes, in principle, two different biological systems. The first is an idealized ecological system of interacting 
species. To represent more realistic ecological systems, the connectivity of the model should be built correspondingly, 
typically, with several levels of preys and predators. 

On the other hand, the model can also describe the system of genotypes present in, or accessible to, a single species 
or population. Within a single species, the number of competing genotypes can be much larger than the number of 
competing species in an ecological niche. Of course, not all of them strive and result finally expressed in the living 
population, and this is precisely the problem we have addressed in this work. 

We have found that the evolution of the system follows complicated orbits in phase space. These orbits drive the 
system from the neighborhood of one of the many equilibria to another, regardless of their stability. Systems with a 
finite population threshold, that may represent more accurately real biological systems, eventually fall into a stable 
equilibrium situation. As time elapses, a variable number of populations become extinct through the interaction with 
the others. In general, more than one species survive, in contrast with the "principle of competitive exclusion" ||] 
(that is known to be of limited validity). The number of surviving ones -those that finally reach equilibrium- is 
characterized by a bell-shaped distribution, whose width and maximum depend on the distribution of interactions. 

Besides competition, a population of genotypes is also subject to changes that arise from random mutations and 
recombination during the reproduction of the organisms [|l3| . The description of such a system would require a 
modification of model, whose behavior cannot be predicted a priori. Mutations can be easily taken into account by 
allowing a new interaction, namely random transitions between the genotypes . The analysis of this system is the 
subject of work in progress. 
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APPENDIX A: 



Theorem. Let {Kij} be an M x M random matrix, whose coefficients are drawn from the same probability 
distribution p{k), with p{k) = for k < 0. Then, the probabihty that the solution to the set of linear equations 

M 

^At,,n, = l, (^ = 1,2,...,M), (Al) 

i=i 

has positive components, > V i = 1, 2, M, is P = 2^^*^. 

Proof. For M = 1, rii — I/kh, which is always positive (P = 1). For M = 2, the system can be explicity solved and, 
in particular, we get 

712 Kii - K21 

The symmetry of this expression with respect to the coefficients Kij makes clear that ni/n2 > with probability 1/2 , 



irrespectively of the form of p{k). Now -since > V i, j- if tii and 712 have the same sign and satisfy Eqs. (Al) 
with M — 2, they must be positive. Therefore, P — 1/2 



For M > 2, we take any pair of equations from (Al) -say the fc-th and the l-th- and rewrite them as 



This can also be put in the form 



Klknk+ Kuril = 1 - I]j^fe^; 



f^'lk^k + k'hUi =1, 



with = Kkk/{^ — J2j=^k I '^kjnj), and analogous expressions for k'^,;, kJ^. and k!^. 

Note that for fixed, arbitrary values of rij {j ^ k,l), the functional form of the primed coefhcents in terms of 
the original ones is the same. The probability distributions for kJ^^., kJ.;, kJj, and k'^ are therefore identical. Hence, 
as a consequence of the previous result for M = 2, the probability that rik/ni is positive equals 1/2. This holds 
irrespectively of the distribution for the primed coefhcients, i.e. irrespe ctiv ely of the values of rij (j 7^ k,l). The 



relative sign of any two components, rik and n/, of the solution to Eqs. (Al) is then statistically independent of the 
values of the other components. 

To insure the positivity of all the components it is sufficient to consider M — 1 ratios rik/ni, for instance, ni/ni 
with I — 2,3,...,Af. According to the above results, the probability that all these ratios are positive is (l/2)*^~^ 



Now -since Kij > V if all the rii have the same sign and satisfy Eqs. (Al), they must be positive. Therefore, 
P= (1/2)^^-1 = 21-^-^. 
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Abramson & Zanette. Figure 1 




FIG. 1. Time evolution of the density of some selected populations, from a system consisting of 20 populations. Both axis 
are logarithmic, to emphasize how some of the populations, after becoming almost extinct, grow again to significative values. 
Inset: the same curves, in a log-linear plot. 
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Abramson & Zanette. Figure 2 
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FIG. 2. Time evolution of the total density. Full line: system without density threshold. Dashed line: the same system 
with density threshold. Note the logarithmic scale in the time axis. 
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Abramson & Zanette. Figure 3 
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FIG. 3. Distribution of the number of survivors, P(s). Each curve corresponds to 2000 roahzations of systems with Ak 
as shown in the legend. Inset: The position of the maximum of the distribution, S, as a function of the width Ak of the 
distribution p{k). 
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Abramson & Zanette. Figure 4 
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FIG. 4. The probability of survival of a single species, Ps, as a function of its initial density, no- The two curves correspond 
to 2000 realizations of a symmetric and an asymmetric system. 



Abramson & Zanette. Figure 5 
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FIG. 5. The probability of survival, Pa as a function of the initial derivative of the density, dn/dt\t=o- The two curves 
correspond to 2000 realizations of a symmetric and an asymmetric system. 
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